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Abstract 

The finite-element approach to lattice field theory is both highly accurate 
(relative errors ~ 1/iV 2 , where N is the number of lattice points) and ex¬ 
actly unitary (in the sense that canonical commutation relations are exactly 
preserved at the lattice sites). In this paper we construct matrix elements 
for the time evolution operator for the anharmonic oscillator, for which the 
continuum Hamiltonian is H = p 2 /2 + \q 2k /2k. Construction of such rna- 
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trix elements does not require solving the implicit equations of motion. Low 
order approximations turn out to be quite accurate. For example, the ma¬ 
trix element of the time evolution operator in the harmonic oscillator ground 
state gives a result for the k = 2 anharmonic oscillator ground state energy 
accurate to better than 1%, while a two-state approximation reduces the er¬ 
ror to less than 0.1%. Accurate wavefunctions are also extracted. Analogous 
results may be obtained in the continuum, but there the computation is more 
difficult, and not generalizable to field theories in more dimensions. 

I. INTRODUCTION 

For over a decade now, the finite-element method has been developed for application to 
quantum systems. (For a review of the program see JIJ.) The essence of the approach is 
to put the Heisenberg equations of motion for the quantum system on a Minkowski space- 
time lattice in such a way as to preserve exactly the canonical commutation relations at 
each lattice site. Doing so corresponds precisely to the classical fini te-element, prescription 
of requiring continuity at the lattice sites while imposing the equations of motion at the 
Gaussian knots, a prescription chosen to minimize numerical error. We have applied this 
technique to examples in quantum mechanics and to quantum field theories in two and 
four space-time dimensions. In particular, recent work has concentrated on Abelian and 
non-Abelian gauge theories n especially on issues of chiral symmetry breaking. 

Because it is the equations of motion that are discretized, a lattice Lagrangian does not 
exist in Minkowski space. This is because the equations of motion are in general nonlocal, 
involving fields at all previous (but not later) times. Similarly, a lattice Hamiltonian does 
not exist, in the sense of an operator from which the equations of motion can be derived. 

However, because the formulation is unitary, a unitary time-evolution operator must 
exist which carries fields from one lattice time to the next. For linear finite elements this 
operator in quantum mechanics has been explicitly constructed [7). Construction of this 
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operator requires solving the equations of motion, which are implicit. Therefore, it is most 
useful, and perhaps surprising, that when matrix elements of the time evolution operator are 
constructed in a harmonic oscillator basis, they do not require the solution of the equations of 
motion j8|. Although these general formulas were derived some years ago, it seems they have 
not been exploited. Our purpose here is to study, in a simple context, the matrix elements 
of the evolution operator, and see how accurately spectral information and wavefunctions 
may be extracted.^ Our goal, of course, is to apply similar techniques in gauge theories, for 
example, to study chiral symmetry breaking in QCD. 

II. REVIEW OF THE FINITE-ELEMENT METHOD 

Let us consider a quantum mechanical system with one degree of freedom governed by 
the continuum Hamiltonian 

H=^ + V(q), (2.1) 

from which follow the Heisenberg equations 

P=-V\q), q=p. (2.2) 

These equations are to be solved subject to the initial condition 

[?(0),p(0)] = i. (2.3) 

It immediately follows from ( |2.2| ) that the same relation holds at any later time 

[q(t),p(t)\=i- (2.4) 

Now suppose we introduce a time lattice by subdividing the interval (0, T) into N subin¬ 
tervals each of length h. On each subinterval (“finite element”) we express the dynamical 
variables as rth degree polynomials 

1 A preliminary version of this work appears in ||. 
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(2.5) 


PW ='}2‘‘k(t/h) k , q(t) = '^2b k (t/h) t , 
k =0 k =0 

where t is a local variable ranging from 0 to h. We determine the 2(r+l) operator coefficients 
a k , b k , as follows: 

1. On the first finite element let 


ao=Po=p(0), b 0 = q 0 = q(0). 


( 2 . 6 ) 


2. Impose the equations of motion (|2.2| ) at r points within the finite element, at cq/r, 
i — 1, 2,..., r, where 0 < aq < «2 < • • • < a r < 1. This then gives 


p(h) gj Pi = a fc , gib) -- 7 1 = y" 


(2.7) 


fc =0 


k=0 


3. Proceed to the next finite element by requiring continuity (but not continuity of deriva¬ 
tives) at the lattice sites, that is, on the second finite element, set 


°o — Pi, bo — qi, 


( 2 , 8 ) 


and again impose the equations of motion at ath, and so on. 

How are the cxj’s determined? By requiring preservation of the canonical commutation 
relations at each lattice site, 

[qi,Pi] = [qo,Po]=h ( 2 - 9 ) 

one finds 


r = 1 (linear finite elements) : 

1 

a = - 
2 

(2.10a) 

: 2 (quadratic finite elements) : 

1 , 1 

a± = ± _ 

2 2y/3 

(2.10b) 

r = 3 (cubic finite elements) : 

1 1 

ai ’ 3_ 2 T 2v / 5’ 012 ~~ 2 

(2.10c) 
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These points are exactly the Gaussian knots, that is, the roots of the rth Legendre polyno¬ 
mial, 


P r (2a - 1) = 0. 


( 2 . 11 ) 


Amazingly, these are precisely the points at which the numerical error is minimized. It is 
known for classical equations that if one uses N rth degree finite elements the relative error 
goes like N~ 2r , while imposing the equations at any other points would give errors like N~ r . 

Let us consider a simple example. The quartic anharmonic oscillator has the continuum 
Hamiltonian 

H=\p 2 + \\q\ ( 2 . 12 ) 

for which the equations of motion are 

q — p, p = —Xq 3 . (2.13) 


If we use the linear (r = 1) finite-element prescription given above, the corresponding discrete 
lattice equations are 


gi ~ Qo _ Pi + Po Pi -Po 
h ~ 2 ’ h 


X 
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(91 + qo) 3 - 


(2.14) 


(Notice the easily remembered mnemonic for linear finite elements: Derivatives are replaced 
by forward differences, while undifferentiated operators are replaced by forward averages.) 
By commuting the first of these equations with pi+po and the second with qi+qo the unitarity 
condition (|2.9| ) follows immediately. These equations are implicit, in the sense that we must 
solve a nonlinear equation to find qi and p\ in terms of q 0 and p 0 . Although such a solution 
can be given, let us make a simple approximation, by expanding the dynamical operators at 
time 1 in powers of h, with operator coefficients at time 0. Those coefficients are determined 
by ( |2.14|) , and a very simple calculation yields 


qi = q 0 + hp 0 - ~h 2 ql + . . . , 

3 

Pi = Po - Xhql - -Xh 2 q 0 poqo + .... (2.15) 
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We can define Fock-space creation and annihilation operators in terms of the initial-time 
operators 


which satisfy 


(cl T (3^) 

q ° = 


V o = 


(a — ad) 
iy/ 2y 


(2.16) 


[a, ad] = 1. 


(2.17) 


Here we have introduced an arbitrary variational parameter 7 , which represents the width 
of the corresponding harmonic oscillator states. The Fock-space states (harmonic oscillator 
states) are created and destroyed by these operators: 


1 \ (ad) n . \ 

\n) = ~7=f l°)> 

Vn! 


(2.18) 


which states are not energy eigenstates of the anharmonic oscillator. We can now take 


matrix elements in these states of the dynamical operators at lattice site 1, using (|2.15|) : 


(l|Pi|0> ~ <l|po|0)(l + i^hX'y 4 - + ...) 

« <l|p o |0)(1 + iujh - ^u; 2 h 2 + ...), (2.19) 

(l|?i|0) ~ (l|go|0)(l + *4 _ 4 2 V + • • •) 

T 4 

« (l|g 0 |0)(l + icoh - 4 ) 2 h 2 + ...), (2.20) 


where we have assumed approximately exponential dependence on the energy difference u. 
Equating the coefficients of the terms through order h 2 constitutes four equations in two 
unknowns. These equations are consistent and yield 

= | a 7 4 = 4 ’ ( 2 - 21 ) 

2 7 ^ 

so the energy difference between the ground state and the first excited state is approximately 

/o x i/3 

uj — ( -A j ~ 1.145A 1//3 (2.22) 
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which is only 5% higher than the exact result E m = 1.08845A 1 / 3 . A similar calculation using 
quadratic finite elements (r = 2) reduces the error to 0.5%. 


Numerical results for a large number of energy differences can also be obtained by taking 
the discrete Fourier transform of the time sequence {(0|g n |l)}. For example, for 1000 finite 
elements, energy differences are computed at the 2-3% level [llj . 


III. THE TIME-EVOLUTION OPERATOR 


Because the canonical commutation relations are preserved at each lattice site, we know 
that there is a unitary time evolution operator that carries dynamical variables forward in 
time: 


Qn+l^UqnUK Vn+\=UVnU\ 


(3.1) 


For the system described by the continuum Hamiltonian (2.1) in the linear finite-element 
scheme, we have found [7] the following formula for U: 


JJ _ e ihpl/4 e ihA(q n ) e ihpl/4: 


(3,2) 


whereQ 

A(x) = ~ 9 ~ l {^x/h 2 )} 2 + V(g~ 1 (4x/h 2 )), (3.3) 

g(x) = j^ x + v '( x ). (3.4) 

The implicit nature of the finite-element prescription is evident in the appearance of the 
inverse of the function g. 

Given the time evolution operator, a lattice Hamiltonian may be defined by U — 
exp(?7i7d). For linear finite elements 7i differs from the continuum Hamiltonian by terms of 
order h 2 . For example, 


2 A misprint occurs in (2.21b) of Q. 
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T r 1 2 2 

V = —m q : 
2 

y = ^ 3: 


74 = —- tan 
mh 


-i 


m/?.\ 


1 2 I 1 2 2 

, -p H —m q 

2/22 


ft = + 7^Ag 3 + /r 2 

ft = ^P 2 + ^Ag 4 + h 2 


X 

12 


pqp + p 


a 2 6 A 2 ■ 

‘24 9 “ 8 qP q 


+ . • •, 
+ 


(3.5a) 

(3.5b) 


(3.5c) 


If one uses quadratic finite elements 7d differs from the continuum Hamiltonian by terms of 
order h 4 , etc. 


IV. MATRIX ELEMENTS OF DYNAMICAL VARIABLES 


Remarkably, it is not necessary to solve the equations of motion to compute matrix 
elements of the dynamical variable. Introduce creation and annihilation operators as in 
( |2.16| ). Then, in terms of harmonic oscillator states (|2.18|) the following formula is easily 
derived || for a general matrix element of q \: 


'y 

(m|gi|n)=- j=(y/m5 ntin -i + A/n5 mjn _i) 


g— i6(m—n ) 


V2 

dzze- 92 ^ 4R2 g\z)H n (g(z)/2R)H m (g(z)/2R), 


(4.1) 


Rs/ tt 2 n+m n\m\ J- c 

where g is given by (p.4|), H n (x) is the nth Hermite polynomial, and we have introduced the 


abbreviations 


R 2 = 


4y 2 1 


+ 


e-“ = i4+ ! 


(4.2) 


ft 4 h 2 7 2 ’ ” Rh 2 ' Rh-i 

For the example of the harmonic oscillator, this formula gives for the ground state-first 
excited state energy difference u> = (2/h) tan _1 (h/2), consistent with (|3.5a|) , while for the 
anharmonic oscillator (|2.12|) if we expand in h we obtain precisely the expansion (|2.20|). 


V. MATRIX ELEMENTS OF THE TIME EVOLUTION OPERATOR 


A similar formula can be derived for the harmonic oscillator matrix elements of the time 


evolution operator. (There is an error in the formula printed in [§.) 




















(m\U\n) = 


1 


1 


2 R \J'K2 n+m n\m\ 


3 — i(n+m+l)9 


x / dzg'{z)H n (g(z)/2R)H m (g(z)/2R)e^ z)+ih3v '^ 2/8 — h29 ^ )2e i6 ^ R] , (5.1) 

J — OO 

which again is expressed in terms of g not g -1 . 

For the harmonic oscillator, where V = q 2 / 2, (|5.1|) gives for the ground-state energy 

(0|f/|0) = e iuJ ° h 


1 _ 1 h 

‘""=ft tan_ 2' 


(5.2) 


which follows from (|3.5a|) . For the general anharmonic oscillator, V = Xq 2k /2k, again we 
expand in powers of h, with the result, for the harmonic oscillator ground state, 


(0|f/|0) = 1 + ihX 1/(k+1) f(a) - yA 2/{k+1) s(a) 
~ 1 + iuj^h — —(Uq h 2 + ..., 


where a = A 7 2A,+2 and 
1 


/(«) = 


s [a = 


Aa l K k ~ l_1 ) 

1 


1 + 


2 «r(fc + 1 / 2 ) 
~k r(l/ 2 ) 


3 Aa 2k ~ ir(fc + l/2) 4a 2 r(2fc + l/2) 


(5.3) 

(5.4a) 

(5.4b) 


lQ a 2 /(k+i) ~ k r(l/ 2 ) ' k 2 T(l/ 2 ) 

This result also derivable from (3.5), using ( 2.16|) , but with considerably more labor. Equat 


ing powers of h in (5.3) gives us two equations, which are to be solved first for the dimen¬ 
sionless number a. Once the number a is determined, the value of Uq is expressed as 


wo = A 1 /( ‘ + 1 >/(a). 


(5.5) 


For a first estimate, we use only the 0(h) equation (|5.5|) and employ the “principle of 
minimum sensitivity” (PMS) [llj] that is, use the stationary value of a, 


f'(a ) = 0 a = 


->k—l 


(2k — 1)!! 


/(«) = 


k + 1 [(2k — l)!!] 1 /^ 41 ) 


4 k 


2 (fc-i)/(fc+i) 


(5.6) 


Specific examples are 


/(«) = < 


0.4293, k = 2, 
0.4639, k = 3, 
0.5230, k = 4, 


(5.7) 
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which are higher than the exact values [T2j of 0.420805, 0.43493, and 0.46450 by about 2%, 
7%, and 13%, respectively. In fact, when we solve ( |5.3| ) for a, that is solve /(a ) 2 = s(a), 
we find complex values, for example, for k = 2 


1 + i 

a = - ± —-= 

2 2 ^ 


/(a) = 0.4178 =F 0.00771. 


(5.8) 


The imaginary part is small, and the real part is only 0.7% low. The corresponding result 
for k = 3 is /(a) = 0.4453 =|= 0.03521, where the real part is now only 2% high. However, for 
k — 4, /(a) = 0.5171 =F 0.07131, and the real part is still 11% high. The failure of (|5.8| ) to be 
real does not indicate any breakdown of unitarity, but only that the one state approximation 
is not exact. 

We do much better by making a two-state approximation, where we must diagonalize 
the 2 x 2 matrix 


% u ^ 

Uoo tv 02 


V 
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U20 U22 

For k = 2 we then find the following relation between uv 0j2 and a = A 7 6 : 

A % 3 

u 0 9 = —-cT 1/3 [12 + 21a =f 2^(8 + 16a + 33a 2 ) 1/2 ], 

16 


(5.9) 


(5.10) 


which, for the — sign, is plotted in Fig. |l|. This graph shows that the ground-state energy 
is very insensitive to the value of a for a broad range of values. The principle of minimum 
sensitivity gives 


u 0 = 0.42124A 1/3 , 


(5.11) 


in spectacular agreement with the exact result, being only 0 . 1 % high, while it gives a value 
for the third state, = 2.992A 1 / 3 , accurate to 1%. (The exact value is 2.959A 1 / 3 [p~3||.) 
Solving for a from the eigenvalues of ( |5.9|) gives even better results: 

uj 0 = A 1/3 (0.42054 ± 2 x 10’%), u ; 2 = A 1 /s (2.9433 ± 0.02201), (5.12) 
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where the ground state energy is now low by 0.06%, the imaginary part being negligible. 
The real part of the energy of the third state is low by only 0.5%. These results for the 
ground state are much better than those resulting from the WKB approximation [13[ . 


The corresponding PMS values for the k — 3 and k — 4 oscillators are similarly im¬ 
pressive: u ; 0 = 0.43913A 1 / 4 (+0.9%) and a ; 0 = 0.47718A 1 / 5 (+2.7%), respectively. Using the 
0{h 2 ) data to compute a gives truly outstanding agreement: 


(Jo = (0.43284 ± 0.00259z)A 1/4 

(5.13a) 

u 2 = (3.4532 ±0.127U)A 1/4 

(5.13b) 

coo = (0.46347 ±0.01463f)A 1/5 

(5.13c) 

= (4.0186 ± 0.308B)A 1/5 

(5.13d) 


where the real parts of the ground-state energies are now low by 0.5% and 0.2%, respectively. 


VI. WAVEFUNCTIONS 


In the process of diagonalizing O we find the corresponding wavefunctions in the two- 
state approximation, that is, the wavefunctions are taken to be linear combinations of n = 0 
and n — 2 harmonic oscillator states of width 7 . The real parts of these wavefunctions are 
plotted in Figs. || and |j. (The imaginary parts amount to only 2% for the ground-state 
wavefunction and 5% for the excited-state wavefunction.) These are normalized to unity at 


the origin to facilitate comparison with [13 


When these are compared with the exact wavefunctions given in [|T3[, we see that the 
approximate ground-state wavefunction is nearly indistinguishable from the exact one, and 
is much better that the WKB wavefunction given there. The excited-state wavefunction 
is quite good, but deviates slightly from the exact wavefunction, and in particular, the 
minimum at x ~ 1.1 should be about 10% deeper. This deviation is not surprising, since the 
exact excited-state wavefunction must contain a substantial mixing with the n = 4 harmonic 
oscillator state. The error in the excited-state wavefunction is also manifested in that fact 
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that in this approximation the two wavefunctions are not quite orthogonal, the magnitude 
of the overlap being 5%. 


VII. CONCLUSIONS 


The simple calculations given here for quantum-mechanical anharmonic oscillators are 
the beginning of a program to develop use of lattice Hamiltonian techniques to explore gauge 
theories in the finite-element context. The numerical results presented in Sec. [Y] and 0 also 


hold in the continuum, by virtue of (|3.5|) , but the calculations are much more tedious without 
the use of the finite-element formula Q5.1|) . It is in two or more space-time dimensions that 
the essential nature of the lattice in such calculations comes into play p].p|.[l4| . The high 
accuracy contrasted with the simplicity of the approach leads us to expect that we can 
extract spectral information, anomalies, and symmetry breaking from an examination of 
the time-evolution operator in gauge theories. 
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FIG. 1. Ground-state energy for the quartic anharmonic oscillator as a function of a = A 7 6 , in 
the second approximation. Here too = A 1 / 3 f(a), f(a) given by ( [5.10 ). 
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